function distance = haversine(lat1, lon1, lat2, lon2)
    % 将经纬度转换为弧度
    lat1 = deg2rad(lat1);
    lon1 = deg2rad(lon1);
    lat2 = deg2rad(lat2);
    lon2 = deg2rad(lon2);
    
    % 应用haversine公式
    dlat = lat2 - lat1;
    dlon = lon2 - lon1;
    a = sin(dlat / 2) * sin(dlat / 2) + cos(lat1) * cos(lat2) * sin(dlon / 2) * sin(dlon / 2);
    c = 2 * atan2(sqrt(a), sqrt(1 - a));
    
    % 地球平均半径，单位为公里
    R = 6371;
    distance = R * c;
end